subroutine symspec(k,nxpol,ndiskdat,nb,nbslider,idphi,nfsample,fgreen,   &
     iqadjust,iqapply,gainx,gainy,phasex,phasey,rejectx,rejecty,         &
     pxdb,pydb,ssz5a,nkhz,ihsym,nzap)

!  k        pointer to the most recent new data
!  nxpol    0/1 to indicate single- or dual-polarization
!  ndiskdat 0/1 to indicate if data from disk
!  nb       0/1 status of noise blanker
!  idphi    Phase correction for Y channel, degrees
!  nfsample sample rate (Hz)
!  fgreen   Frequency of green marker in I/Q calibrate mode (-48.0 to +48.0 kHz)
!  iqadjust 0/1 to indicate whether IQ adjustment is active
!  iqapply  0/1 to indicate whether to apply I/Q calibration
!  pxdb     power in x channel (0-60 dB)
!  pydb     power in y channel (0-60 dB)
!  ssz5a    polarized spectrum, for waterfall display
!  nkhz     integer kHz portion of center frequency, e.g., 125 for 144.125
!  ihsym    index number of this half-symbol (1-322)
!  nzap     number of samples zero'ed by noise blanker

  parameter (NSMAX=60*96000)          !Total sample intervals per minute
  parameter (NFFT=32768)              !Length of FFTs
  real*8 ts,hsym
  real*8 fcenter
  common/datcom/dd(4,5760000),ss(4,322,NFFT),savg(4,NFFT),fcenter,nutc,junk
  real*4 ssz5a(NFFT),w(NFFT)
  complex z,zfac
  complex zsumx,zsumy
  complex cx(NFFT),cy(NFFT)
  data rms/999.0/,k0/99999999/,nadjx/0/,nadjy/0/
  save

  if(k.lt.NFFT) then
     ihsym=0
     go to 999             !Wait for enough samples to start
  endif
  if(k0.eq.99999999) then
     pi=4.0*atan(1.0)
     do i=1,NFFT
        w(i)=(sin(i*pi/NFFT))**2
     enddo
  endif

  hsym=2048.d0*96000.d0/11025.d0      !Samples per JT65 half-symbol
  if(nfsample.eq.95238)   hsym=2048.d0*95238.1d0/11025.d0
  npts=NFFT                           !Samples used in each half-symbol FFT

  if(k.lt.k0) then
     ts=1.d0 - hsym
     savg=0.
     ihsym=0
  endif
  k0=k

  ihsym=ihsym+1
  ja=ts+hsym                          !Index of first sample
  jb=ja+npts-1                        !Last sample

  ts=ts+hsym
  px=0.
  py=0.
  i=0
  fac=0.0002
  dphi=idphi/57.2957795
  zfac=fac*cmplx(cos(dphi),sin(dphi))
  nzap=0
  sigmas=1.5*(10.0**(0.01*nbslider)) + 0.7
  peaklimit=sigmas*max(10.0,rms)
  do j=ja,jb                          !Copy data into cx, cy
     x1=dd(1,j)
     x2=dd(2,j)
     if(nxpol.ne.0) then
        x3=dd(3,j)
        x4=dd(4,j)
     else
        x3=0.
        x4=0.
     endif
     if(nb.ne.0) then
        peak=max(x1,x2,x3,x4,-x1,-x2,-x3,-x4)
        if(peak.gt.peaklimit) then
           x1=0.
           x2=0.
           x3=0.
           x4=0.
           dd(1:4,j)=0.
           nzap=nzap+1
        endif
     endif
     i=i+1
     cx(i)=fac*cmplx(x1,x2)
     cy(i)=zfac*cmplx(x3,x4)          !NB: cy includes dphi correction
     px=px + x1*x1 + x2*x2            !Compute power
     py=py + x3*x3 + x4*x4
  enddo
  if(nzap/178.lt.50 .and. (ndiskdat.eq.0 .or. ihsym.lt.280)) then
     nsum=jb-ja+1-nzap
     if(nsum.le.0) nsum=1
     rmsx=sqrt(0.5*px/nsum)
     rmsy=sqrt(0.5*py/nsum)
     rms=rmsx
     if(nxpol.ne.0) rms=sqrt((px+py)/(4.0*nsum))
  endif
  pxdb=0.
  pydb=0.
  if(rmsx.gt.1.0) pxdb=20.0*log10(rmsx)
  if(rmsy.gt.1.0) pydb=20.0*log10(rmsy)
  if(pxdb.gt.60.0) pxdb=60.0
  if(pydb.gt.60.0) pydb=60.0
!  if(mod(ihsym,10).eq.0) write(*,3009) ihsym,nzap,rmsx,rmsy,rms,pxdb,pydb
!3009 format(i4,i8,5f10.2)

  cx=w*cx                             !Apply window
  if(nxpol.ne.0) cy=w*cy

  call four2a(cx,NFFT,1,1,1)          !Compute FFTs
  if(iqadjust.eq.0) nadjx=0
  if(iqadjust.ne.0 .and. nadjx.lt.50) call iqcal(nadjx,cx,gainx,phasex,     &
                                                 zsumx,ipkx,rejectx0)
  if(iqapply.ne.0) call iqfix(cx,NFFT,gainx,phasex)

  if(nxpol.ne.0) then
     call four2a(cy,NFFT,1,1,1)
     if(iqadjust.eq.0) nadjy=0
     if(iqadjust.ne.0 .and. nadjy.lt.50) call iqcal(nadjy,cy,gainy,phasey,  &
                                                 zsumy,ipky,rejecty)
     if(iqapply.ne.0) call iqfix(cy,NFFT,gainy,phasey)
  endif

  if(nadjx.eq.50) then
!     freqx=ipkx*96.0/nfft - 48.0
!     freqy=ipky*96.0/nfft - 48.0
     rejectx=rejectx0
!     write(*,3011) freqx,gainx,phasex,rejectx,                   &
!                   freqy,gainy,phasey,rejecty
!3011 format('X  Freq:',f9.3,'   Gain:',f8.4,'   Phase:',f9.5,'   Image:',f6.1/&
!            'Y  Freq:',f9.3,'   Gain:',f8.4,'   Phase:',f9.5,'   Image:',f6.1)
     nadjx=9999
     nadjy=9999
  endif

  n=ihsym
  do i=1,NFFT
     sx=real(cx(i))**2 + aimag(cx(i))**2  
     ss(1,n,i)=sx                    ! Pol = 0
     savg(1,i)=savg(1,i) + sx
     
     if(nxpol.ne.0) then
        z=cx(i) + cy(i)
        s45=0.5*(real(z)**2 + aimag(z)**2)
        ss(2,n,i)=s45                   ! Pol = 45
        savg(2,i)=savg(2,i) + s45

        sy=real(cy(i))**2 + aimag(cy(i))**2
        ss(3,n,i)=sy                    ! Pol = 90
        savg(3,i)=savg(3,i) + sy
        
        z=cx(i) - cy(i)
        s135=0.5*(real(z)**2 + aimag(z)**2)
        ss(4,n,i)=s135                  ! Pol = 135
        savg(4,i)=savg(4,i) + s135

        z=cx(i)*conjg(cy(i))
        q=sx - sy
        u=2.0*real(z)
        ssz5a(i)=0.707*sqrt(q*q + u*u)    !Spectrum of linear polarization
! Leif's formula:
!     ssz5a(i)=0.5*(sx+sy) + (real(z)**2 + aimag(z)**2 - sx*sy)/(sx+sy)
     else
        ssz5a(i)=sx
     endif
  enddo
  if(ihsym.eq.278) then
     if(iqadjust.ne.0 .and. ipkx.ne.0 .and. ipky.ne.0) then
        rejectx=10.0*log10(savg(1,1+nfft-ipkx)/savg(1,1+ipkx))
        rejecty=10.0*log10(savg(3,1+nfft-ipky)/savg(3,1+ipky))
!        write(*,3012) rejectx,rejecty
!3012    format('After I/Q balance   X Image:',f6.1,'dB   Y Image:',f6.1,' dB')
     endif
!     df=96.0/NFFT
!     do i=1,NFFT
!        freq=df*(i-1) - 48.0
!        write(71,3101) freq,(10.0*log10(savg(j,i)),j=1,4)
!3101    format(f12.6,4f12.3)
!     enddo
!     call flush(71)
  endif

  nkhz=nint(1000.d0*(fcenter-int(fcenter)))
  if(fcenter.eq.0.d0) nkhz=125

999 return
end subroutine symspec
